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PREFACE 


The  Elliptic  Restricted  Three-Body  Problem  (ER3BP)  continues  to  be 
the  focus  of  an  increasing  level  of  attention  in  the  scientific 
community.  Equilibrium  points  defined  within  a  localized  system  of 
gravitational  and  centrifugal  forces  in  a  rotating  three-body  system 
were  the  subject  of  considerable  research  during  the  past  200  years. 
More  recently,  spacecraft  orbits  near  these  equilibrium  points  have 
been  of  great  interest;  in  particular,  the  ER3BP  defined  to  include  a 
spacecraft  in  orbit  within  the  gravitational  attractions  of  the  Sun  and 
the  Earth  has  proven  to  be  of  great  value  in  the  study  of  the 
solar-terrestrial  environment.  The  addition  of  solar  radiation 
pressure  to  the  gravitational  and  centrifugal  iorce  structure  axters 
the  locations  and  conditions  for  linear  stability  of  the  equilibrium 
points  in  the  Sun-Earth  ER3BP.  (These  equilibrium  points  are  also 
sometimes  referred  to  as  libration  points  or  Lagrange  points.)  Solar 
radiation  pressure  also  consequently  affects  the  spacecraft  orbits 
computed  near  these  points. 

This  work  explores  the  generation  of  bounded  orbits  in  the  ER3BP, 
and  some  of  the  results  of  adding  the  solar  radiation  pressure  force  to 
the  ER3BP  are  presented.  The  focus  is  limited  to  the  equilibrium 
points  defined  in  the  Sun-Earth  ER3BP,  although  extensions  to  other 
systems  may  be  straightforward.  One  equilibrium  point  in  this 
three-body  system  is  located  between  the  Sun  and  the  Earth;  orbits  near 
this  point  are  the  primary  focus  for  the  study  of  solar  effects  on  the 
terrestrial  environment  and  are  therefore  the  orbits  considered  here. 
This  effort  is  supported  by  the  Frank  J.  Seiler  Research  Laboratory  and 
has  been  conducted  as  doctoral  research  under  the  direction  o' 
Professor  K.C.  Howell,  School  of  Aeronautics  and  Astronautics,  Purdue 
University,  West  Lafayette,  Indiana. 
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CHAPTER  1:  INTRODUCTION 


With  the  expansion  of  space  exploration  programs  worldwide, 
interest  has  increased  in  the  design  of  innovative,  complex,  and  yet 
low-cost  spacecraft  trajectories  that  meet  demanding  mission 
requirements.  In  most  of  the  missions  flown  in  the  last  few  decades, 
the  spacecraft  spent  the  majority  of  the  flight  time  in  a  force 
environment  dominated  by  a  single  gravitational  field.  For  Llv* 
preliminary  mission  analysis  in  these  cases,  additional  attracting 
bodies  and  other  forces  could  be  modeled,  when  required,  as  perturbing 
influences.  Analysis  of  some  recently  proposed  and  more  adventurous 
missions,  such  as  those  involving  libration  point  orbits,  will  require 
dynamic  models  of  higher  complexity,  since  at  least  two  gravitational 
fields  are  of  nearly  equal  influence  on  the  spacecraft  throughout  the 
majority  of  the  mission.  Thus,  trajectories  determined  for  a  system 
consisting  of  numerous  gravitational  forces  have  been  of  particular 
theoretical  and  practical  interest  in  recent  years. 

One  type  of  many-body  problem,  motion  within  a  three-body  system 
of  forces,  has  a  wide  range  of  applications.  The  general  problem  of 
three  bodies  assumes  that  each  body  has  finite  mass  and  that  the  motion 
is  a  result  of  mutual  gravitational  attraction.  When  the  mass  of  one 
of  the  three  bodies  is  assumed  to  be  sufficiently  small  (infinitesimal) 
so  that  it  does  not  affect  the  motion  of  the  other  two  bodies 
(primaries)  in  the  system,  the  "restricted  three-body  problem"  results. 
The  primaries  can  be  further  assumed  to  be  moving  in  known  elliptic  or 
circular  orbits  about  their  common  center  of  mass.  Therefore,  the 
elliptic  restricted  three-body  problem,  where  the  primaries  are  assumed 
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to  be  in  known  elliptic  orbits,  may  be  considered  a  reasonably 
approximate  model  for  a  spacecraft  moving  within  the  gravitational 
fields  of  the  Sun  and  the  Earth,  for  instance. 


A.  Definition  of  the  Problem 

In  the  formulation  of  the  restricted  three-body  problem,  one  mass 
is  defined  as  infinitesimal  relative  to  the  remaining  two  masses 
(primaries).  The  primaries,  unaffected  by  the  infinitesimal  mass,  move 
under  their  mutual  gravitational  attractions.  In  the  elliptic 
restricted  three-body  problem  (ER3BP),  the  primaries  are  assumed  to 
move  on  elliptic  paths.  If  the  eccentricity  of  the  primaries’  orbit  is 
assumed  to  be  zero,  the  circular  restricted  three-body  problem  (CR3BP) 
results.  Even  for  known  primary  motion,  however,  a  general, 
closed-form  solution  for  motion  of  the  third  body  of  infinitesimal  mass 
does  not  exist.  In  the  restricted  three-body  problem  (ER3BP  or  CR3BP), 
five  equilibrium  (libration)  solutions  can  be  found.  These  equilibrium 
points  are  particular  solutions  of  the  equations  of  motion  governing 
the  path  of  the  infinitesimal  mass  moving  within  the  gravitational 
fields  of  the  primaries. 

The  equilibrium  points  are  defined  relative  to  a  coordinate  system 
rotating  with  the  primaries.  At  these  locations,  the  forces  on  the 
spacecraft  are  in  equilibrium.  These  forces  include  the  gravitational 
forces  from  the  massive  bodies  and  the  centrifugal  force  associated 
with  the  rotation  of  the  system.  (The  addition  of  solar  radiation 
pressure  to  the  force  model  changes  the  locations  of  the  five  Lagrange 
points,  although  they  can  still  be  defined,  and  these  solar  radiation 
effects  will  be  discussed  in  the  next  chapter.)  The  five  libration 
points  are  located  in  the  plane  of  primary  rotation  and  are  depicted  in 
Figure  1-1  for  the  CR3BP.  Three  of  the  libration  points  are  on  the 
line  between  the  two  massive  bodies,  and  one  of  these  collinear  points 
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is  interior  to  the  primaries.  The  last  two  points  are  at  the  vertices 
of  two  equilateral  triangles  in  the  plane  of  primary  rotation.  The 
triangles  have  a  common  base  that  is  the  line  between  the  primary 
masses. 

For  the  CR3BP,  the  five  libration  points  are  stationary  relative 
to  the  rotating  reference  frame.  The  approximate  relative  distances 
between  the  primaries  and  the  five  libration  points  are  also  given  in 
Figure  1-1.  If  the  problem  is  generalized  to  the  ER3BP,  the  libration 
points  pulsate  as  the  distance  between  the  primaries  varies  with  time. 
In  both  the  circular  and  the  elliptic  restricted  problems, 
two-dimensional  and  three-dimensional  trajectories,  both  periodic  and 
quasi-periodic  paths,  can  be  computed  in  the  vicinity  of  these 
libration  points. 

B.  Previous  Contributors 

The  formulation  and  solution  of  the  equations  of  motion  governing 
two  or  more  heavenly  bodies  moving  under  the  influence  of  their  mutual 
gravitational  fields  have  inspired  many  researchers  through  the 
centuries.  The  difficulty  involved  in  finding  any  solution  to  the 
multi-body  problem  has  been  fortunate  for  the  development  of  the 
mathematics,  science,  and  engineering  methods  used  in  the 
investigations.  A  few  of  the  many  significant  contributions  are 
mentioned  briefly  here.  The  problem  history  is  described  in  terms  of 
the  many-body  problem,  the  three-body  problem,  and  then  libration  point 
orbits  that  include  both  periodic  and  quasi-periodic  trajectories. 

1.  The  Many-Body  Problem 

Even  the  simplified  two-body  problem  successfully  investigated  by 
Johanness  Kepler  was  not,  in  fact,  simple  to  solve.  Kepler’s  laws  of 
planetary  motion,  published  between  1609  and  1619,  did  not  account  fci 
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the  masses  of  the  two  planets  and  assumed  that  the  planets  move  on 
elliptical  paths  about  the  Sun.  The  findings  of  Kepler  were  a  major 
development  in  that  they  supported  the  heliocentric  theory,  which 
documented  evidence  suggests  was  first  postulated  by  Aristarchus  of 
Samos  in  Ancient  Greece  and  then  rediscovered  by  Nicholas  Copernicus 
(1473-1543). 

Sir  Isaac  Newton  (1642-1727)  was  first  to  formulate  the  many-body 
problem  as  three  or  more  massive  particles  moving  under  their  mutual 
gravitational  attractions.  Other  researchers,  including  Kepler,  had 
proposed  the  concept  of  the  inverse  square  law  of  gravitation,  but 
Newton’s  many  scientific  discoveries  included  his  recognition  that 
the  inverse  square  law  of  gravitation  was  the  fundamental  law  of 
celestial  mechanics.  However,  even  with  the  planets  modeled  as  point 
masses,  finding  any  solution  to  the  problem  involving  three  or  more 
bodies  is  nontrivial. 

Leonhard  Euler  (1707-1783)  showed  that  ten  integrals  of  the  motion 
(ten  constants  of  integration)  exist  in  the  many-body  problem.  Six  of 
these  scalar  constants  result  because  total  linear  momentum  is 
conserved.  Three  of  the  scalar  constants  of  integration  result  from 
the  conservation  of  total  angular  momentum.  The  tenth  constant  is  due 
to  the  conservation  of  total  mechanical  energy  of  the  system  of  bodies. 
Euler  was  the  first  to  use  mathematical  methods,  rather  than  the 
accepted  geometrical  approaches  of  his  time,  to  solve  problems  in 
dynamics,  and  he  is  therefore  widely  considered  to  be  the  father  of 
analytical  dynamics.  Many  years  later,  Henri  Poincare  (1854-1912) 
significantly  concluded,  along  with  others,  that  no  additional 
integrals  exist  in  the  many-body  problem  if  three  conditions  hold:  (1) 
the  integrals  of  the  motion  are  equations  involving  only  integral  or 
algebraic  functions  of  the  coordinates  and  velocities  of  the  bodies  in 
question,  (2)  they  are  valid  for  all  values  of  the  masses,  and  (3)  they 
satisfy  the  equations  of  motion.  Poincare  had  extended  and  generalized 
the  work  of  H.  Bruns  whose  findings  were  published  in  1887  and 
described  specific  necessary  conditions  for  the  existence  of  integrals 
of  the  motion  for  the  general  problem  of  three  bodies. 
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The  set  of  differential  equations  governing  motion  in  the 
many-body  problem  includes  a  second-order  vector  differential  equation 
for  each  body.  When  n  bodies  are  defined,  3n  second-order  scalar 
differential  equations  result.  In  order  to  completely  specify  the 
motion  of  these  n  particles,  6n  constants  of  integration  are  required. 
Since  Euler  showed  that  10  integrals  exist  and  Poincare  proved  that 
only  the  10  (with  certain  restrictions  on  form)  can  exist,  a  three-body 
problem  that  requires  18  constants  of  the  motion  would  thus  defy 
general  closed-form  solution.  By  eliminating  time  as  the  independent 
variable  and  by  using  an  "elimination  of  the  nodes"  method  originally 
developed  by  Carl  Gustav  Jacobi  (1804-1851),  the  many-body  problem  may 
be  reduced  to  order  6(n-2).  The  general  three-body  problem  would  still 
have  six  unknown  constants  of  the  motion.  It  is  intriguing  that  the 
general  three-body  problem  looks  deceptively  similar  to  the  general 
two-body  problem,  yet  the  general  problem  of  three  bodies  is  of  a  much 
higher  degree  of  complexity  and  does  not  have  a  general  closed-form 
solution. 


2.  The  Three-Body  Problem 

The  three-body  problem  has  long  been  a  rich  research  area  in 
celestial  mechanics.  Newton  studied  and  found  some  approximate 
solutions  to  the  general  problem  of  three  bodies.  In  1765,  Euler 
solved  the  problem  of  three  massive  bodies  that  maintain  a  constant 
ratio  of  their  relative  distances  from  each  other  and  move  on  a  fixed 
straight  line.  Later,  Joseph-Louis  Lagrange  presented  his  particular 
solutions  of  the  general  three-body  problem  in  1772.  Lagrange,  a 
student  and  friend  of  Euler,  found  five  equilibrium  (libration)  points 
in  the  general  three-body  problem.  These  five  libration  points  are 
thus  also  referred  to  as  Lagrange  points.  Pierre-Simon  de  Laplace 
repeated  Lagrange’s  analysis  in  his  works  published  between  1799  and 
1825.  Joseph  Liouville  (1809-1882)  studied  the  Sun-Earth-Moon 
three-body  problem  in  particular.  Liouville  investigated  the 
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instability  of  the  three-body  configuration  and  the  effects  of 
disturbing  forces  on  the  Moon’s  motion.  He  proved  that  the  collineo 
libration  (or  Lagrange)  points  are  unstable  points  of  equilibrium.  The 
two  triangular  points  have  been  found  to  be  linearly  stable  for  certain 
ratios  of  the  primary  masses  and  for  certain  values  of  primary  orbit 
eccentricity.  George  W.  Hill  also  studied  the  motion  of  the  Moon  in 
the  Sun-Earth-Moon  system  and  published  his  theory  of  lunar  motion  ir. 
1378. 

The  "restricted"  problem  of  three  bodies  has  been  of  great 
research  interest  in  the  last  100  years.  Poincare  was  apparently  the 
first  researcher  to  coin  the  term  "restricted"  to  refer  to  the  problem 
in  which  the  infinitesimal  mass  is  assumed  to  have  no  eifect  on  the 
movement  of  the  finite  masses.  The  primary  masses  may  be  assumed  to  be 
in  either  circular  or  elliptic  orbits,  and  the  motion  of  the 
infintesimal  mass  in  the  vicinity  of  a  Lagrange  point  can  then  be 
studied. 


3.  Libration  Point  Orbits 

In  this  section,  a  brief  historical  background  of  libration  point 
orbital  research  is  presented.  Then,  one  goal  of  the  continuing 
scientific  investigations  concerning  libration  point  orbits  in  the 
Sun-Earth  system  is  discussed.  The  third  part  of  this  section 
describes  several  planned  missions  in  the  Sun-Earth  system.  The  final 
two  subjects  discussed  are  periodic  (halo)  and  quasi-periodic 
(Lissajous)  orbits,  respectively. 


a.  Historical  Background 


Libration  point  orbits  have  been  the  subject  of  considerable 
research  interest.  Poincare  completed  significant  work  on  bounded 
periodic  orbits  that  were  two-dimensional  (orbits  of  his  first  and 
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second  kind)  in  the  CR3BP  and  the  ER3BP  (respectively)  and  that  were 

three-dimensional  (orbits  of  his  third  kind)  in  the  CR3BP.  He  looked 

at  families  of  periodic  orbits  indexed  by  the  ratio  of  the  mass  of  the 

smaller  primary  to  the  sum  of  the  primary  masses.  Works  by  Sir  George 

Howard  Darwin  (son  of  Charles  Darwin)  in  1899  and  1911  investigated  the 

circular  restricted  three-body  problem.  His  studies  included  the 

stipulations  that  the  infinitesimal  mass  remained  in  the  plane  of 

motion  of  the  primaries  and  that  one  primary  was  ten  times  as  massive 
[1  2] 

as  the  other.  ’  Da "win  numerically  calculated  a  number  of  periodic 
orbits  in  the  planar  CR3BP.  (Of  course,  numerical  integration 
computations  at  that  time  were  highly  labor-intensive.)  Karl  Frithiof 
Sundman  obtained  a  general  solution  using  convergent  series  for  the 
general  problem  of  three  bodies  in  1912,  but  the  slowness  of 

[3  4) 

convergence  made  it  of  little  practical  value.  ’  Henry  Crozier 
Plummer  published  works  in  1903  and  1914  concerning  approximate 
analytic  solutions  to  the  planar  CR3BP.  Forrest  Ray  Moulton 

described  the  problem  in  his  1914  book  and  analyzed  a  great  variety  of 
periodic  orbits  and  approximate  analytic  solutions  in  the  planar  CR3BP 
in  1920.  [7,8) 

Several  additional  scientific  discoveries  and  research 

accomplishments  in  this  century  helped  to  spur  interest  in  libration 

point  orbits.  Until  the  early  twentieth  century,  investigations 

concerning  the  five  Lagrangian  points  had  been  considered  interesting 

14] 

but  purely  academic.  It  had  seemed  highly  unlikely  that  such 

unusual  planetary  formations  could  exist  in  nature.  Then,  in  1906,  a 
small  planet  in  the  Sun-Jupiter  three-body  system  was  found  to  be 
moving  in  close  proximity  to  one  of  the  triangular  libration  points. 
This  small  planet  is  named  Achilles  and  is  one  of  the  Trojan  asteroids 
that  now  have  been  found  to  number  approximately  25.  Before  the  1950s, 
solutions  were  also  constrained  to  hand  calculations  and  analytic 
approximations.  A  1966  book  by  Victor  Szebehely  thoroughly  summarized 
the  evolution  of  efforts  to  understand  the  three-body  problem  to  that 

[9] 

date.  Until  then,  most  of  the  work  of  generating  solutions  had  been 
limited  to  two-dimensional  orbits  in  the  plane  of  motion  of  the 
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prim?"ies.  Within  the  past  30  years,  high  speed,  electronic  computing 
devices  helped  make  more  accurate  numerical  solutions  possible  and 
enabled  the  study  of  three-dimensional  trajectories. 

More  recently,  the  discovery  of  other  objects  orbiting  in  the 
vicinity  of  libration  points  has  also  helped  to  increase  interest. 
Analysis  of  data  from  the  Voyager  1  spacecraft  found  that  two  of 
Saturn’s  moons  complete  the  Saturn-Dione  A-Dione  B  equilateral  triangle 
solution  for  that  three-body  system.  The  two  moons  are  the  same 
distance  from  Saturn,  and  the  much  smaller  Dione  B  slowly  oscillates 

o 

near  a  point  that  remains  some  60  ahead  of  Dione  A  as  the  moons  orbit 

around  Saturn.  Some  meteoric  dust  particles  are  also  suspected  to  be 

orbiting  near  the  triangular  libration  points  in  the  Earth-Moon 
(4  9] 

system.  ’  In  the  Sun-Earth  system,  the  collinear  libration  point  on 
the  dark  side  of  the  Earth  may  be  the  location  of  meteoric  dust 
particles  that  can  be  seen  to  reflect  light  from  the  setting 
Sun.  t4’7’9*  xhe  discoveries  of,  and  speculations  about,  these  natural 
formations  located  in  the  vicinity  of  Lagrange  points,  and  the 
increasing  capabilties  of  electronic  computers,  have  facilitated 
interest  in  orbits  for  artificial  satellites  near  these  points. 

The  circular  or  elliptic  motion  of  the  primaries,  as  modeled  in  the 
CR3BP  or  ER3BP,  may  be  considered  a  good  orbital  approximation.  The 
orbit  of  a  planet  about  the  Sun  or  the  orbit  of  the  Moon  about  the 
Earth  is  approximately  Keplerian,  and  the  spacecraft  can  easily  be 
considered  an  infinitesimal  mass  in  comparison.  The  study  of  libration 
point  orbits  in  such  systems  is  gaining  in  significance  because  the 
positioning  of  spacecraft  n=>ar  Earth-Moon  or  Sun-Earth  libration  points 
can  be  of  value  to  meet  a  number  of  objectives.  Earth-Moon  equilateral 
libration  points  have  been  studied  for  potential  stationing  of  space 
colonies  and  transportation  relay  stations.  The  translunar  libration 
point  in  the  Earth-Moon  system  has  also  been  proposed  as  the  location 
of  a  communications  relay  satellite  to  support  a  lunar  exploration  base 
on  the  far  side  of  the  Moon.  The  Sun-Earth  interior  libration  point  is 
a  site  of  great  interest  for  upstream  investigation  of  solar  effects  on 
the  Earth.  One  mission,  the  International  Sun-Earth  Explorer-3, 
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successfully  conducted  such  scientific  research  while  in  a  libration 
point  orbit  for  nearly  4  years.  Scientific  interest,  then,  continues 
to  motivate  Sun-Earth  libration  point  trajectory  studies;  in  part, 
because  solar  activity  has  some  remarkable — and  damaging — influence  on 
the  Earth’s  atmosphere. 


b.  One  Objective  of  Sun-Earth  Libration  Point  Orbital  Missions 

One  important  goal  of  some  Sun-Earth  libration  point  orbital 
research  is  to  study  solar  effects  on  the  terrestrial  environment.  The 
ceaseless  solar  wind  contains  a  constant  stream  of  particles  stripped 
from  the  Sun;  however,  it  is  the  unexpected  solar  storms  that  do  the 
greatest  amount  of  harm  on  Earth.  The  unpredictable  solar  eruptions 
release  sudden  surges  of  particles  into  the  solar  wind,  and  the 
particles  reach  the  Earth  some  4  days  later  (they  travel  at 
approximately  1.5  million  kilometers  per  hour).1101  The  surge  in 
particles  can  cause  a  geomagnetic  storm  that  upsets  delicate 
electronics  and  that  can  damage  electrical  supply  systems.  Changes  in 
the  Earth’s  magnetic  field  (caused  by  the  geomagnetic  storm)  can  induce 
quasi-direct  currents  in  high-voltage  alternating  current  transmission 
lines  and  can  thus  damage  transformers. 1101  In  March  1989,  a  surge  of 
particles  hit  the  areas  near  the  North  Pole  and  caused  a  blackout  in 
Canada’s  province  of  Quebec  and  the  shutdown  of  a  nuclear  power  plant 
in  New  Jersey.  Two  hundred  electrical  utilities  were  in  some  measure 
affected  by  this  solar  storm. 1101 

Significantly,  there  are  additional  major  terrestrial  changes  as  a 
result  of  the  unpredictable  solar  storms.  Increased  solar  radiation 
also  causes  the  Earth’s  outer  atmosphere  to  expand.  The  outer 
atmospheric  expansion,  in  turn,  increases  the  drag  on  some  satellites, 
thus  degrading  their  attitude  control  and  changing  their  orbital 
tracks.  In  1989,  the  United  States  Air  Force  Space  Command  temporarily 
lost  track  of  approximately  1,400  space  objects  for  up  to  several  weeks 
as  increased  drag  due  to  solar-induced  terrestrial  atmospheric 
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The 


expansion  changed  their  previously-predicted  paths. 1101 
solar-induced  atmospheric  expansion  near  the  Earth  has  also  caused  the 
Soviet  Salyut  7  space  station  to  plunge  to  Earth  in  1991 ,  some  3  years 
early.  The  Soviets  had  boosted  it  to  a  higher  orbit  from  where  it 
should  have  continued  orbiting  until  approximately  the  year  1994; 
however,  the  increased  atmospheric  drag  changed  its  orbital  path  and 
caused  it  to  begin  an  uncontrolled  tumble. 

Unfortunately,  these  solar  storms  are  not  well  predicted. 
Ground-based  solar  storm  warning  systems  are  inadequate;  in  fact, 
currently  only  14%  of  the  storms  are  preceeded  by  an  accurate  warning, 
and  the  existing  warning  system  provides  a  78%  false  alarm  rate.1101 
One  of  the  several  objectives  of  future  libration  point  orbit  research 
in  the  Sun-Earth  system  could  then  be  the  study  of  solar  storms  and  the 
potential  use  of  storm  warning  spacecraft.  Positioning  of  a  warning 
spacecraft  near  the  Sun-Earth  interior  Lagrange  point  will  allow  about 
an  hour’s  notice  of  an  incoming  solar  storm. 


c.  Some  Missions  Planned  in  the  Sun-Earth  System 

Future  spacecraft  missions  will  increase  scientific  knowledge  of 

the  solar-terrestrial  environment  and  are  planned  as  collaborative 

efforts  of  several  agencies.  The  National  Aeronautics  and  Space 

Administration  (NASA),  the  European  Space  Agency  (ESA),  and  the 

Institute  of  Space  and  Astronautical  Science  (ISAS)  in  Japan  are  the 

major  participants,  with  some  minor  involvement  of  the  United  States’ 

Department  of  Defense  (DoD),  in  this  solar-geospace  research.  1111 

Several  planned  missions  i;e  part  of  the  joint  International 

Solar-Terrestrial  Physics  (UIP)  program.  A  mission  designated  as  WIND 

is  scheduled  for  launch  on  December  22,  1992;  following  a  2-year 

double-lunar  swingby  trajectory,  WIND  will  be  injected  into  a  libration 

112] 

point  orbit  between  the  Earth  and  the  Sun  in  January  1995.  The 

Solar  and  Heliospheric  Observatory  (SOHO)  Mission  is  currently 

1 121 

scheduled  for  launch  on  March  14,  1995.  SOHO  will  be  placed  into 
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orbit  near  the  1 ibration  point  between  the  Earth  and  the  Sun,  and  its 
related  scientific  research  activities  are  planned  to  continue  for  6 
years. 

Other  missions,  coordinated  with  and  scientifically  complementing 
the  vehicles  in  the  libration  point  orbits,  will  be  launched  in  the 
next  decade,  as  well.  A  set  of  four  spacecraft  called  CLUSTER  will  be 
located  in  polar  Ear/  .  bits  at  the  same  time  that  SOHO  is  in  its 
libration  po5nt  K\  .  CLUSTER  will  provide  complementary 

solar-terrestriai  . ;e •*  •,.<  .roni  inside  the  Earth’s  magnetosphere. 1111 
A  NASA-controlled  mis  ion,  POLAR,  is  planned  for  injection  into  a 
northern  Earth  polar  orb5 t  in  order  to  take  measurements  of  the 
magnetosphere.  GEOT--'’  an  fSAS-control led  mission,  is  planned  to 
primarily  study  the  g*  magnetic  tail.  A  joint  NASA/DoD  mission,  the 
Combined  Release  and  Radiation  Effects  Satellite  (CRRES),  is  planned  to 
be  in  an  Earth  orbit  in  order  to  investigate  the  near-Earth  plasma 
sheet  and  the  outer  Van  Allen  belt. 1111  Numerous  specific  missions  are 
also  being  planned  for  launch  into  libration  point  orbits  after  the 
turn  of  the  century.  Scientific  interest  in  libration  point  orbits  of 
various  sizes  and  shapes,  tailored  to  specific  missions,  will  thus 
continue  to  increase. 


d.  Halo  0-  bits 

Three-dimensional,  periodic  "halo"  orbits  in  the  vicinity  of  the 
collinear  libration  points  have  been  studied  since  the  late  1960s. 
Early  work  concerning  these  orbits  was  motivated  by  studies  related  to 
exploring  the  far  side  of  the  Moon.  These  studies  were  completed  in 
support  of  the  planned  Apollo  1C  lunar  exploration  mission  that  was 
later  canceled.  Robert  Farquhar  coined  the  term  "halo"  to  describe  a 
three-dimensional,  periodic  orbit  near  a  libration  point  on  the  far 

[13] 

side  of  the  Moon  in  the  Earth-Moon  system.  These  orbits  were 

designed  to  be  large  enough  so  that  the  spacecraft  would  be  constantly 
in  view  of  the  Earth  and  would  thus  appear  as  a  halo  around  the  Moon. 
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A  communications  station  in  this  type  of  orbit  could  maintain  constant 
contact  between  the  Earth  and  a  lunar  experimentation  station  on  the 
far  side  of  the  Moon. J14) 

Interest  in  these  periodic,  three-dimensional  orbits  has 
continued.  Since  the  late  1960s,  a  number  of  researchers  have  produced 
approximate  analytic,  as  well  as  numerically  integrated,  solutions  that 
represent  three-dimensional,  periodic  halo  orbits  in  the  restricted 
three-body  problem;  some  have  also  conducted  stability  studies  on  many 

(15-22] 

families  of  such  orbits.  In  the  CR3BP,  precisely  periodic  halo 
paths  can  be  constructed,  and  the  orbital  symmetry  of  the  constructed 
trajectory  can  be  fully  exploited  to  facilitate  the  numerical 
computations.  However,  in  the  ER3BP.  orbital  construction  is  more 
difficult,  and  halo-type  orbits  that  are  net  perfectly  periodic  may 
result. 

An  orthographic  projection  of  approximately  one  revolution  of  a 
halo-type  orbit  in  the  Sun-Earth  ER3BP  is  depicted  in  Figure  1-2.  This 
view  of  the  orbit  is  only  one  of  the  three  possible  orthographic 
depictions  of  a  three-dimensional  orbit,  and  it  is  included  here  for 
illustrative  purposes — the  actual  orbits  used  in  this  research  will  be 
fuliy  described  in  the  following  chapter.  The  orbit  is  generateo  near 
the  libratic..  point  located  between  the  Sur  and  the  Earth.  This  planar 
projection  is  perpendicular  to  the  line  between  the  primaries;  thus,  it 
is  the  view  of  the  orbit  as  seen  when  looking  towara  the  Sun  from  the 
Earth.  Note  that  the  z  coordinate  represents  out-cf-plan*1  .notion. 
Since  the  distance  between  the  Sun  and  the  Earth  is  app'  >.  r  |,®ly 

3 

1.5x10  kilometers,  the  z  excursion  is  relatively  small  but 
significant.  Note  also  that  the  maximum  y  excursion  for  this  halo-type 
orbit  is  nearly  four  times  as  large  as  its  maximum  z  excursion. 

A  recently  successful  mission  in  which  the  spacecraft  moved  on  a 
halo  trajectory  about  a  libration  point  between  the  Earth  and  the  Sun 
has  helped  spur  continuing  interest  in  libration  point  orbits.  On 
August  12,  1978,  a  scientific  spacecraft  named  the  International 
Sun-Earth  Explorer-3  (ISEE-3)  was  launched  toward  a  halo  orbit  near  the 
interior  Sun-Earth  libration  point.  It  remained  in  orbit  for  nearly  4 
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z  excursion  from  the  libration  point  in  km 


Ar  indicates  direction  of  motion 


Figure  1-2.  An  Orthographic  Depiction  of  a  Halo-Type  Orbit. 
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years  to  measure  solar  wind  characteristics  and  for  the  study  of  other 
solar-induced  phenomena  such  as  solar  flares  upstream  of  the  Earth. 
Measurements  associated  with  these  processes  could  be  made  at  the 
interior  Lagrange  point  approximately  1  hour  before  any  disturbance 
reached  the  terrestrial  environment.  Concurrently,  the  ISEE-1  and 
ISEE-2  missions  measured  similar  characteristics  from  Earth  orbit  in 
order  to  complete  a  comparative  study.  Future  missions  planned  as  part 
of  the  upcoming  ISTP  program  are  being  designed  to  continue  the 
investigation  of  the  solar-terrestrial  environment.1111 


e.  Lissajous  Orbits 

Quasi-periodic  orbits  near  libration  points  are  also  currently  of 
great  research  interest.  One  orthographic  view  of  several  revolutions 
of  a  Lissajous  orbit  constructed  near  the  interior  libration  point  in 
the  Sun-Earth  ER3BP  is  depicted  in  Figure  1-3.  This  is  the  view  of  the 
orbit  as  seen  when  looking  toward  the  Sun  from  the  Earth.  This  single 
view  of  the  three-dimensional  Lissajous  orbit  is  included  for 
illustrative  purposes  and  for  a  comparison  with  the  halo-type  orbit 
depicted  in  Figure  1-2.  Note  that  the  z  excursions  of  the  two  orbits 
are  comparable;  yet,  the  maximum  y  excursion  of  the  Lissajous  path  is 
only  approximately  one-fourth  that  of  the  halo-type  orbit. 

The  variations  in  size  and  shape  that  a  quasi-periodic  orbit  can 
exhibit  may  add  valuable  flexibility  for  mission  planning.  This  type 
of  bounded,  three-dimensional  libration  point  trajectory  is  called  a 
Lissajous  orbit  since  specific  planar  projections  of  these 
quasi-periodic  trajectories  may  look  like  a  special  type  of  "Lissajous" 
curve.  Physicist  Jules  Antoine  Lissajous  (1822-1880)  investigated 
curves  that  were  generated  by  compounding  simple  harmonic  motions  at 
right  angles,  and  he  delivered  a  paper  on  this  subject  to  the  Paris 
Academy  of  Sciences  in  1857.  Nathaniel  Bowditch  of  Salem, 
Massachusetts,  had  conducted  some  similar  work  in  1815.  Lissajous 
curves  have  a  wide  variety  of  shapes  that  depend  on  the  frequency, 
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indicates  direction  of  motion 


Figure  1-3.  An  Orthographic  Depiction  of  a  Lissajous  Orbit. 
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[23  24] 

phase,  and  amplitude  of  the  orthogonal  components  of  the  motion.  ’ 

When  the  in-plane  and  the  (orthogonal)  out-of-plane  frequencies  of  the 
linearized  motion  are  nearly  (but  not)  equal,  the  resulting  path  is 

w* 

typically  called  a  Lissajous  trajectory. 

A  method  to  generate  approximations  for  this  type  of 

quasi-periodic  orbital  path  was  developed  analytically  by  Farquhar  and 
[151 

Kamel  in  1973.  They  derived  a  third-order  approximate  analytic 

solution  for  a  translunar  libration  point  orbit  in  the  Earth-Moon  ER3BP 

that  also  included  solar  gravity  perturbations.  In  1975,  Richardson 

and  Cary  then  developed  a  fourth-order  analytic  Lissajous  approximation 

in  the  Sun-Ear th+Moon  barycenter  system.  The  notation  "Earth+Moon" 

indicates  that  the  Earth  and  the  Moon  are  treated  as  one  body  with  mass 

center  at  the  Earth-Moon  barycenter.  In  consideration  of  the  lunar 

perturbation,  Farquhar  has  shown  that  the  accuracy  of  solutions  in  the 

Sun-Earth  CR3BP  can  be  enhanced  if  the  collinear  libration  points  are 

defined  along  the  line  between  the  Sun  and  the  center  of  mass  of  the 

[17] 

Earth  and  the  Moon.  Since  1975,  a  few  researchers  have  considered 

methods  to  numerically  generate  Lissajous  trajectories,  but  the  lack  of 
periodicity  of  a  Lissajous  path  complicates  numerical  construction  of 
bounded  trajectories.  Howell  and  Pernicka  have  developed  a  numerical 
technique  for  determination  of  three-dimensional,  bounded  Lissajous 
trajectories  of  arbitrary  size  and  duration.112’26”301 


C.  Overview 

This  next  chapter  considers  the  research  results  and  briefly 
discusses  the  following  topics:  coordinate  systems,  equations  of 
motion,  locations  of  the  five  Lagrange  points,  linearization  of  the 
equations  in  order  to  compute  the  state  transition  matrix,  and 
construction  of  bounded  orbits  near  the  equilibrium  points.  Follow-on 
studies  (not  presented  here)  can  then  be  used  to  conduct  tracking  and 
control  scheme  simulations,  incorporating  the  solar  radiation  pressure 
force,  for  a  spacecraft  in  such  an  orbit. 
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CHAPTER  2:  RESULTS 


In  this  chapter,  the  elliptic  restricted  three-body  problem  and 
the  associated  coordinate  systems  are  reviewed;  the  equations  of  motion 
for  an  infinitesimal  mass  moving  in  the  gravity  fields  of  two  massive 
bodies  are  then  derived.  Locations  of  the  libration  points  are 
discussed,  and  then  the  computation  of  the  state  transition  matrix  and 
the  construction  of  bounded  nominal  orbits  near  the  collinear  Lagrange 
points  are  summarized. 


A.  Elliptic  Restricted  Three-Body  Problem 

The  elliptic  restricted  three-body  problem  is  a  simplification  of 
the  general  problem  of  three  bodies.  In  the  general  three-body 
problem,  each  of  the  three  bodies  is  assumed  to  be  a  particle  of  finite 
mass  and,  thus,  exerts  an  influence  on  the  motion  of  each  of  the  other 
bodies.  Neither  the  general  nor  the  restricted  problem  of  three  bodies 
has  a  general  closed-form  solution.  However,  when  problem 
simplifications  are  made,  particular  solutions  can  be  determined.  If 
the  mass  of  one  of  the  bodies  is  restricted  to  be  infinitesimal,  such 
that  it  does  not  affect  the  motion  of  the  other  two  massive  bodies 
(primaries),  the  restricted  three-body  model  results.  The  primaries 
are  assumed  to  be  in  known  elliptic  (ER3BP)  or  circular  (CR3BP)  orbits 
about  their  common  mass  center  (barycenter).  The  problem  can  then  be 
completely  described  by  a  single  second-order  vector  differential 
equation  with  variables  appropriately  defined  for  a  specified 
coordinate  frame. 
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B.  Coordinate  Systems 


The  two  standard  coordinate  systems  used  in  the  analysis  of  this 

problem  have  a  common  origin  at  the  center  of  mass  (barycenter)  of  the 

primaries.  Primaries  with  masses  m  and  m  such  that  m  £  m  are 

assumed  here,  although  this  distinction  is  arbitrary.  The 

infinitesimal  mass  is  denoted  as  m  .  These  masses  (m  ,m  ,m  ) 

3  12  3 

correspond  to  particles  situated  at  points  P  ,  P2>  and  P^, 

respectively.  The  barycenter  is  denoted  as  "B,“  and  the  resulting 

arrangement  is  shown  in  Figure  2-1.  The  rotating  coordinate  system  is 

defined  as  x  y  z  ,  and  the  inertial  system  is  identified  as  x  y  z  . 

Note  that  both  coordinate  systems  are  right-handed,  and  the  x  and  y 

axes  for  both  systems  are  in  the  plane  of  motion  of  the  primaries.  The 

Xj  axis  is,  of  course,  assumed  to  be  oriented  in  some  fixed  direction; 

in  this  specific  formulation  of  the  problem,  it  is  assumed  to  be 

parallel  to  a  vector  defined  with  a  base  point  at  the  Sun  and  directed 

toward  periapsis  of  the  Earth’s  orbit.  The  rotating  xr  axis  is  defined 

along  the  line  that  joins  the  primaries  and  is  directed  from  the  larger 

toward  the  smaller  primary.  The  z  axes  are  coincident  and  are  directed 

parallel  to  the  primary  system  angular  momentum  vector.  The  yR  axis 

completes  the  right-handed  x  y  z  system. 

R  R  R 


C.  Equations  of  Hotion 

Newtonian  mechanics  are  used  to  formulate  the  equations  of  motion 
for  m^  (the  spacecraft)  relative  to  B  as  observed  in  the  inertial 
reference  frame.  The  sum  of  the  forces  on  m^  resulting  from  both  the 
gravity  fields  of  masses  (the  Sun)  and  m^  (the  Earth-Moon 
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Figure  2-1.  Coordinate  Systems  With  Barycenter  Origin. 
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barycenter)  and  from  the  solar  radiation  pressure  can  be  used  to 
produce  the  following  second-order  vector  differential  equation: 


// 

P 


m 

=  -  G  (— i-)  d  -  G  ( 
d3 


m 


2_ 

3 


r 


)  r  + 


(2-1) 


The  overbar  denotes  a  vector,  and  primes  indicate  differentiation 

with  respect  to  dimensional  time.  All  quantities  are  dimensional,  as 

appropriate,  and  the  quantity  "G"  is  the  universal  gravitiational 

constant.  The  scalars  "d"  and  "r"  in  equation  (2-1)  denote  the 

magnitudes  of  the  vectors  d  and  r,  respectively,  as  depicted  in 

Figure  2-1.  The  dimensionless  scalar  "k"  is  the  solar  reflectivity 

constant,  and  "S"  is  the  solar  radiation  pressure  constant.  The 

formulation  of  the  solar  radiation  force  model  and  the  values  for  the 

[  31 1 

solar  radiation  constants  are  derived  from  previous  work  by  Bell. 

The  numerical  values  for  these  constants  are  selected  from 
characteristic  data  for  ISEE-3;  however,  they  should  also  approximately 
model  a  vehicle  such  as  SOHO  in  the  International  Solar-Terrestrial 
Physics  Program: 


S  = 


A  S  r2 
o  o 

c  m  ’ 


(2-2) 


where  k  =  1.2561,  Sq=  1352.098  kg/sec3,  r  =  1 . 4959787x10s  km,  c  =  the 

g 

speed  of  light  =  2.998x10  m/sec,  A  =  surface  area  of  the  spacecraft’s 

2 

sun-facing  side  =  3.55  m  ,  and  m3  =  mass  of  the  spaceciaft  =  435  kg. 

The  constant  Sq  is  the  solar  light  flux  measured  at  one  astronomical 

unit  (A.U.),  and  r  is  the  nominal  distance  associated  with  the 

131 1 

measurement  of  the  solar  flux  S  .  The  value  of  A,  also  termed  the 

effective  cross-sectional  area  of  the  spacecraft,  is  considered  to  be 
constant  in  this  work.  The  constant  k  is  a  material  parameter 
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dependent  on  the  absorptivity  of  the  spacecraft  surface  and  is 

[31 ) 

generally  confined  to  a  range  of  0  i  k  i  2.0. 

The  position  vector  p  is  written  in  rotating  components  as 


p  =  x  x  +  y  y  +  z  z  (2-3) 

K  r  J  JR  R 

where  x  ,y  ,z  are  unit  vectors.  The  velocity  and  the  acceleration  of 
r  ■'r  r  j 

the  spacecraft  (particle  P3  with  mass  m3)  relative  to  the  barycenter  B 
as  observed  in  the  inertial  reference  frame  can  then  be  derived  by 
using  the  basic  kinematic  expressions 

P  -  p  =  p  +  w  x  p,  (2-4) 


and 


"Jt 

P 


I— R 

W 


x  p 


(2-5) 


where  the  superscript  "i"  denotes  the  inertial  reference  frame,  the 

j 

superscript  "r"  indicates  the  rotating  reference  frame,  p  denotes 

—  ^  I 

the  time  derivative  of  p  in  equation  (2-4)  taken  with  respect  to  the 

I*~R 

rotating  frame,  and  w  is  the  angular  velocity  of  the  rotating 

reference  frame  relative  to  the  inertial  frame.  The  symbol  "x" 
indicates  a  vector  product. 

Assuming  an  elliptic  primary  orbit  and  a  known  time  of  periapsis 
passage  for  the  primary  orbit,  two-body  orbital  mechanics  can  be  used 
to  determine  the  position  of  P2  in  its  orbit  at  any  given  time.  In  the 
standard  two-body  problem,  the  angular  displacement  of  d2  from 

periapsis  is  defined  by  the  angle  0  as  depicted  in  Figure  2-1. 

Clearly,  the  angle  0  in  this  particular  formulation  is  also  the  angular 
displacement  of  the  rotating  reference  frame  relative  to  the  fixed 

inertial  system. 
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Hence,  the  angular  velocity  of  the  rotating  frame  relative  to  the 
inertial  direction  is  written 

ISR  =  0'z.  (2-6) 

Using  equations  (2-3)  <nd  (2-6)  in  equation  (2-4),  an  expression  for 
the  velocity  of  P^  can  be  obtained: 

p'  =  p'1  =  (x'  -  0'y)  x  +  (y'  +  0'x)  y  +  z'z  .  (2-7) 

K  K  3  R  3  R 

Then,  using  equations  (2-6)  and  (2-7)  in  equation  (2-5),  the  following 
kinematic  expression  for  p"  can  be  derived: 

p"  =  (x"-0"y-20' y' -0' 2x)x  +(y"+0"x+20' x' -O' 2y )y  +z"z  .  (2-8) 

R  RR 

Notice  that  equation  (2-8)  includes  terms  with  the  angular 
velocity,  0'(t),  and  angular  acceleration,  0"(t),  of  the  rotating 
reference  frame  with  respect  to  the  inertial  frame.  By  using  known 
elliptic  orbital  elements,  readily  computed  expressions  for  0'  and  0" 
can  be  found.  For  instance,  the  magnitude  of  the  specific  angular 
momentum  vector  can  be  computed  by  using  a  vector  equation  that 
includes  0' : 

h  =  hz  =  Rxv=Rx  x  (R'x  +  Ro'y  )  =  rV£  (2-9) 

R  R  R 

where  the  vector  v  is  the  orbital  velocity  of  the  smaller  primary,  P  , 
in  its  elliptic  orbit,  as  observed  in  the  fixed  frame.  For  the 
elliptic  two-body  orbit,  the  varying  distance  R  between  the  primaries 
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is  given  by 


R  =  R(t)  =  a  (1  -  e  cos(E)), 


(2-10) 


where  "e"  is  the  eccentricity,  "E"  is  the  eccentric  anomaly,  and  "a"  is 
the  semi-major  axis  of  the  primary  orbit. 

For  any  orbit,  the  semi-latus  rectum,  p,  is  expressed 


P 


h2/  G  ( m  +  m  ) . 
1  2 


(2-11) 


For  any  conic  orbit  that  is  not  a  parabolic  path,  p  can  also  be  written 


p  =  a(l  -  e)2  =  a(l-e)(l+e). 


(2-12) 


Using  equations  (2-11)  and  (2-12),  an  expression  for  h  can  be  found: 


h  =  Va(l-e) ( 1+e)  VG(m  +  m  ) .  (2-13) 

1  2 

Hence,  using  equations  (2-9),  (2-10),  and  (2-13),  the  angular  rate,  6', 
can  be  computed  as 


,  V'a(l-e)  (1+e)  v'Gtm  +m  ) 
h  _  _ l  2 

R2  a2  (1  -  e  cos(E))2 


(2-14) 


Equation  (2-14)  includes,  however,  the  angle  "E"  that  varies  with 
time.  Therefore,  the  expression  for  the  angular  acceleration,  0",  will 
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necessarily  include  E'  ,  for  which  a.  relationship  in  terms  of  known 
orbital  elements  must  be  found.  An  equation  that  relates  the  eccentric 
anomaly  and  the  time  since  periapsis  passage  is  given  as 


t  -  t 


■(/ 


a3/G(m  + 


m2) 


(E  -  e  sin  (E) ) 


(2-15) 


where  t  is  the  time  of  periapsis  passage.  The  quantity  / G(mj+  m  )/&3 
is  denoted  as  the  mean  motion,  n  ;  that  is, 


mean 


n 

mean 


=  /G(mi+  mz)/a3  , 


(2-16) 


and  equation  (2-15)  can  be  rewritten  as 


n  (t  -  t  )  *  (E  -  e  sin  (E) ) .  (2-17) 

mean  p 


Equations  (2-15)  or  (2-17)  can  be  used  to  derive  a  mathematical 
relationship  involving  the  time  derivative  of  E  and  other  known  orbital 
elements.  Equation  (2-17)  will  later  prove  to  be  quite  useful  in 
determining  the  angle  E  if  both  the  time,  t,  and  the  time  of  periapsis 
passage,  t  ,  are  known.  Differentiating  equation  (2-15)  with  respect 

p 

to  time  produces 


/ G(m  +m  )/a 

E'  =  - - - — - - .  (2-18) 

a(l-e  cos(E) ) 
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Taking  the  time  derivative  of  equation  (2-14)  and  using  equation 
(2-18),  the  angular  acceleration  of  the  rotating  reference  frame 
relative  to  the  inertial  frame  can  be  computed  as 


„  -2G(m  +m  Ji/d-e)  ( l+e)  e  sin(E) 

9  =  - - 

a J(1  -  e  cos(E) )4 


(2-19) 


Then  equations  (2-1),  (2-2),  (2-8),  (2-14),  (2-15),  and  (2-19) 

comprise  the  governing  set  of  equations  for  motion  of  the  infinitesimal 

mass  within  the  force  environment  of  the  Earth+Moon  (m  )  and  the  Sun 

2 

(m^),  including  the  solar  radiation  pressure  force.  Note  that  equation 

(2-15)  is  listed  here.  Equation  (2-15),  or  its  mathematical  equivalent 

(2-17),  known  as  Kepler’s  equation,  is  used  to  compute  the  angle  E  when 

the  time  since  periapsis  passage  is  known.  Of  a  number  of  different 

procedures,  Newton’s  method  for  finding  the  zeros  of  a  function  proves 

adequate  for  this  purpose  in  this  work;  that  is,  given  t  and  time  t, 

p 

the  correct  value  of  the  angle  E  must  be  computed.  Newton’s  method  is 
iterative  and,  in  this  work,  finds  the  zero  of  a  function  defined  as 


f (E)  =  n  (t  -  t  )  -  E  +  e  sin(E). 

mean  p 


(2-20) 


The  algorithm  reduces  to  finding  iterative  solutions  based  on  the 

..  (26) 
equation 


(t  -  t  )  -  E 

p 


+  e  sin(E  ) 

n-l  n-1 


E  =  E 


mean 


(2-21) 


The  iteration  process  using  equation  (2-21)  is  defined  as  complete 
when  the  condition  |f(E  s  10  12  is  satisfied.  (The  symbol  "  |  |" , 
of  course,  indicates  that  the  absolute  value  is  to  be  taken. )  Thus,  by 
using  equation  (2-21),  the  value  of  the  eccentric  anomaly,  E,  at  any 
time,  t,  relative  to  the  time  of  periapsis  passage,  t  ,  can  be 

p 

calculated  and  then  used  in  equations  (2-14)  and  (2-19)  to  determine  Q' 
and  Q". 

Three  scalar  equations  of  motion  for  can  be  derived  using  the 
dimensional  equations  (2-1),  (2-2),  (2-8),  (2-14),  and  (2-19);  however, 
for  convenience,  the  following  scaling  factors  are  typically 
introduced: 

(1)  The  sum  of  the  masses  of  the  primaries  equals  one 

mass  unit,  (m  +  m  =1  unit  of  mass) 

1  2 

(2)  The  mean  distance  between  the  primaries  equals  one 
unit  of  distance,  (a  =  1  unit  of  distance) 

(3)  The  universal  gravitational  constant  is  equal  to  one 
unit  by  proper  choice  of  characteristic  time. 

(characteristic  time  =  1/n  ;  thus  G  =  1  unit) 

mean 

The  dimensional  equations  of  motion  can  be  simplified  and  scaled 
by  introducing  the  characteristic  quantities  defined  above  and  by 
introducing  the  nondimensional  mass  ratio  fi,  "psuedo-potential"  U,  and 
the  scaled  solar  radiation  constant  s: 


and 


m 

2 


m  +  m 
1  2 


(2-22) 


U  = 

d 


1 

2 


02(x2  + 


a  .  k  S 

y  )  -  — 


(2-23) 
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where  the  dot  denotes  the  derivative  with  respect  to  nondimensional 
time.  The  scaled  solar  radiation  constant,  s,  is  derived  from  the 
dimensional  radiation  constant  denoted  as  S  in  equation  (2-2)  and,  by 
using  the  characteristic  quantities  described  above,  its  value  is 
calculated  as  s  =  6.206597029461384xl0~6  nondimensional  units.  Then, 
the  vector  magnitudes,  "d"  and  "r,"  are  written  in  terms  of  scaled 
quantities  as: 


d  =  [<x  +  (x  R}2  +  y2  +  z2]1/2,  (2-24) 

r  =  ((x  -  R  +  v  R)2  +  y2  +  z2]1/2.  (2-25) 


The  three  scalar  second-order  differential  equations  that  result 
can  be  written  in  terms  of  characteristic  quantities  as 


x  -  2  9  y  =  — —  +  6  y  =  U  +  0  y, 

ox  X 


y  +  2  9  x  =  -T—  -  0  x  =  V  -  9  x, 
oy  y 


au 

az 


(2-26) 


(2-27) 


(2-28) 


where 
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(1-ji)  (x+R(i)  /i(x-R+pR) 


ks(x+Rfi) 


(2-29) 


U 

X 


02x 


U 

y 


(l-ji)y 


liy  ksy 

—  +  —  +  0  y  , 
_3  j3 


(2-30) 


( 1 -p )  z  pz  ksz 


(2-31) 


The  scaled  equations  for  the  angular  velocity  and  angular 
acceleration  of  the  rotating  frame  relative  to  the  inertial  frame  also 
simplify: 


,  V(l-e) (1+e) 

G  =  ~  =  - ,  (2-32) 

r2  (1  -  e  cos(E))2 


-2vTT-eTTl+eT  e  sin(E) 

0  =  - .  (2-33) 

(1  -  e  ccs(E))4 


If  the  primaries  are  assumed  to  be  moving  in  a  circular  orbit, 
then  0  =  0,  R  =  0=  1,  and  equations  (2-26),  (2-27),  (2-28)  reduce  to 
three  scalar  equations  in  the  simplified  form: 
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(2-34) 


where 


x  -  2  y  = 


=  U  , 

X 


y  +  2  x  = 


=  U  , 
y 


(2-35) 


=  U  , 

Z 


(2-36) 


(l-fi)(x+/i)  /i(x-  1+m)  ks(x+fi) 

U  - - - - -  +  - - -  +  x  , 


(2-37) 


(l-|i)y  My  ksy 
U  - - -  -  +  y  , 

y  ,3  3  ,3 

y  d  r  d 


(2-38) 


(l-fji)z  mz  ksz 


U  - - 

z  ,3 


3  A 3 

r  d 


(2-39) 


The  system  of  equations  (2-34),  (2-35),  and  (2-36)  admits  an 

integral  of  the  motion  that  was  originally  found  by  Jacobi: 


J  =  2  U  -  x2  -  y2  -  z2  =  constant. 


(2-40) 
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In  addition  to  offering  special  insight  into  the  CR3BP,  Jacobi’s 
constant  can  be  used  as  a  check  during  numerical  integration  of  the 
equations  of  motion;  however,  it  is  not  a  constant  in  the  ER3BP. 

The  scalar  equations  (2-26),  (2-27),  and  (2-28)  corresponding  to 
the  elliptic  restricted  problem  or  equations  (2-34),  (2-35),  and  (2-36) 
derived  for  the  circular  restricted  problem  can  be  used  to  locate  the 
five  libration  points  in  the  rotating  reference  frame 


D.  Locations  of  the  Lagrangian  Points 

By  using  scalar  equations  (2-34),  (2-35),  and  (2-36)  for  motion  in 
the  CR3BP,  the  locations  of  the  stationary  equilibrium  points  can  be 
determined.  Equations  (2-26),  (2-27),  and  (2-28)  can  be  used  to 
determine  rat' os  of  distances  that  are  constant  in  the  ER3BP.  These 
ratios  are  related  to  the  locations  of  libration  points  that  have,  been 
defined  in  the  ER3BP  and  that  "pulsate"  with  respect  to  the  rotating 
reference  frame  as  the  distance  between  the  primaries  varies  with  time. 


1.  The  CR3BP 

In  the  CR3BP,  the  five  libration  points  are  equilibrium  points  and 
are  stationary  with  respect  to  the  rotating  coordinate  frame,  that  is, 
they  are  locations  at  which  the  forces  on  the  third  body  sum  to  zero. 
The  arrangement  of  points  and  the  corresponding  nondimensional 
distances  are  depicted  in  Figure  2-2.  Note  that  three  of  the  libration 
points  (Lj,  L2>  L3)  are  collinear  with  the  primaries;  one  r.ol linear 
point  (L  )  is  interior  to  the  primaries.  The  retaining  two  points  (L 
and  Lg)  are  located  at  the  vertices  of  two  equilateral  triangles  that 
are  in  the  plane  of  primary  rotation  and  that  have  a  common  base 
between  the  primaries. 
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In  the  CR3BP,  the  libration  points  are  stationary  in  the  rotating 
coordinate  frame.  Stationary  points  are  defined  as  points  at  which  the 
relative  velocity  and  acceleration  are  zero,  such  that 

x  =  y  =  z  =  x  =  y  =  z  =  0.  (2-41) 

Using  equations  (2-41)  in  equations  (2-34)  through  (2-39),  the  useful 
conditions  U  =  U  =  U  =  0  are  found.  The  three  collinear  libration 

x  y  z 

points  can  be  readily  located  by  further  noting  that  y  =  z  =  0  for  the 
points  located  on  the  rotating  xr  axis. 

Then,  the  remaining  equation  that  allows  conditions  on  the 
coordinate  x  to  be  found  is  (2-37).  For  the  collinear  libration 
points,  equation  (2-37)  then  reduces  to 


U  =  0 

x 


(1— p)(x+p) 

,3 


p(x-(l-p))  ks(x+p) 


...  + 


+  X. 


(2-42) 


By  substituting  the  appropriate  coordinates  into  equation  (2-42)  for 
each  of  the  three  collinear  Lagrange  points,  the  locations  of  the  three 
points  can  then  be  individually  determined. 

However,  before  deriving  the  equations  that  are  used  to  locate 
each  of  the  libration  points,  the  values  of  the  constants  used  to  solve 
equation  (2-42)  are  noted  here.  The  value  of  p  used  in  this  work  is 
3.0404234945077  x  10  6  (consistent  with  the  Sun/Ear th+Moon  system). 
The  characteristic  distance  equals  the  semi-major  axis  of  the  primary 
orbit;  that  is,  a  =  1.4959787066  x  108  kilometers.  When  the  solar 
radiation  pressure  force  is  included  in  the  model,  the  values  of  the 
nondimensional  constants  include  the  values  of  k  =  1.2561  and 
s  =  6.206597029461384  x  10  6  scaled  units.  Also,  note  that  when  the 
constants  "k"  or  "s"  are  defined  to  be  zero,  equation  (2-42)  reduces  to 
the  equation  derived  for  the  CR3BP  without  the  solar  radiation  force; 
it  may  be  useful  and  interesting  to  compare  the  positions  of  the 
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collinear  points  with  and  without  solar  radiation  pressure  forces 
included  in  the  model.  First,  however,  it  is  appropriate  to  briefly 
discuss  the  selection  of  the  parameter  values  used  in  this  solar 
radiation  pressure  model. 

The  values  chosen  for  computation  of  the  solar  radiation  pressure 
force  deserve  special  discussion  hei  e  because  they  have  a  minor,  yet 
important,  impact  on  several  of  the  numerical  results.  The  computation 
of  the  solar  radiation  pressure  constant  S  is  described  in  equation 
(2-2).  (The  scaled  solar  radiation  pressure  constant  is  denoted  as  s. ) 
The  values  of  the  parameters  used  in  that  equation  are  selected  from 

[31  ] 

research  by  Bell  using  post-flight  data  for  the  ISEE-3  mission  near 

Lagrange  point  Certainly,  many  factors  can  affect  the  accuracy  of 

this  solar  radiation  pressure  approximation.  The  orientation  of  the 

spacecraft  may  alter  the  area  of  its  Sun-facing  side  (the  variable  A  in 

equation  (2-2)),  and  the  position  of  the  spacecraft  at  a  libration 

point  other  than  may  require  modification  of  at  least  the  variable 

r  in  equation  (2-2). 
o 

The  value  1.2561  chosen  for  the  solar  reflectivity  constant  k 
conforms  with  other  research157,59’801  concerning  ISEE-3;  yet,  it  too 
is  clearly  an  approximation.  The  values  used  in  these  works  for  ISEE-3 
include  estimated  values  for  k  of  1.20  and  1.30;  also,  using  actual 
flight  data  for  ISEE-3,  the  constant  k  was  at  times  treated  as  a 
solved-for  parameter.  While  its  value  varied  widely  during  flight, 
some  of  the  solved-for  values  at  various  epochs  were  1.258,  1.261,  and 
1 .260. 159,801  (It  may  be  interesting  to  note  that  the  reflectivity 
properties  of  any  spacecraft  surface  may  change  with  age,  that  is,  with 
time  on  station,  due  in  part  to  possible  particulate  impacts  and  solar 
radiation.)  Early  work  by  Bell  used  a  k  value  of  1.2561;  this  is  the 
nominal  value  of  k  used  to  generate  the  reference  trajectories  for  the 
orbit  determination  error  analysis  and  station-keeping  studies  in  this 
work.  More  recent  work  by  Bell  includes  a  k  value  of  1.20.  The 
results  of  using  this  smaller  value  of  k  when  determining  the  locations 
of  the  Lagrange  points  will  be  included  in  the  next  few  paragraphs. 
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The  positive  constant  y  is  the  distance  between  the  smaller 
primary  mass  and  the  libration  point  as  depicted  in  Figure  2-2.  For 
libration  point  L^,  x  =  L  =  1  -  p  -  y  (see  Figure  2-2),  and  the 
following  quintic  is  derived  from  equation  (2-42): 


y^  -  (3-p)y*+  (3-2p)y^  -  (p+ks)^  +  2py^  p  =  0.  (2-43) 


The  nondimensional  position  of  L  relative  to  F  in  the  CR3BP  can 

1  2 

be  determined  by  finding  the  single  real  root  of  equation  (2-43),  which 

for  this  problem  with  k  =  1.2651  has  a  value  y  =  .01001184941667781. 

When  solar  radiation  pressure  is  not  included,  the  position  of  Lj  is 

determined  by  the  value  y  =  .01001097733285880.  Hence,  the  position 

of  Lagrange  point  Lt  moves  approximately  130.46188  kilometers  closer  to 

the  Sun  (m  at  P  )  when  solar  radiation  pressure  forces  are  included  in 

[31 ) 

the  circular  model.  Bell  computes  a  value  of  124.60  kilometers  for 
the  change  in  the  position  of  Lj  if  the  solar  radiation  pressure  force 
with  k  =  1.20  is  included. 

For  libration  point  L,  x=L  =  1  -  p  +  y  (see  Figure  2-2),  and 
the  following  quintic  is  derived  from  equation  (2-42): 


^  +  (3-p)y*+  (3-2p)y^  -  (p-ks)y^  -  2py  -  p  =  0.  (2-44) 

6  6  o  2  2 


When  solar  radiation  pressure  forces  with  k  =  1.2561  are  included  in 

the  model,  the  location  of  Lagrange  point  L  is  defined  by  the 

nondimensional  value  y^  =  .01007738021243130.  When  solar  radiation 

pressure  is  not  included,  the  position  of  Lagrange  point  Lg  is 

calculated  as  y  =  .01007824054648482.  Hence,  the  location  of  L  in 
2  2 

the  CR3BP  moves  approximately  128.70414  kilometers  closer  to  the  Sun 
when  solar  radiation  pressure  forces  with  k  =  1.2561  are  included  in 
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[31  ] 

the  model.  Bell  computes  a  value  of  122.90  kilometers  for  the 

*  i* 1  i«..  ¥ 

change  in  the  position  of  L if  the  solar  radiation  pressure  force  with 
k  =  1.20  is  incorporated.  (Note  that  solar  radiation  pressure  forces 
may  not  actually  influence  a  spacecraft  in  an  orbit  that  remains  on  the 
dark  side  of  the  Earth,  that  is,  in  the  Earth’s  shadow.  In  fact,  the 
size  and  shape  of  the  orbit  will  be  reflected  in  the  magnitude  of  the 
solar  radiation  force  on  the  spacecraft  near  L2> ) 

For  libration  point  L3>  x  =  L  =  -  ji  -  y3  (see  Figure  2-2),  and 
the  following  quintic  is  derived  from  (2-42): 


^+(2+p)?%-(l+2/i)^-(l-p-ks)^-(2-2fi-2ks)r3-(l-p-ks)  _  0.  (2-45) 


When  solar  radiation  pressure  with  k  =  1.2561  is  included  in  the  force 
model,  the  location  of  Lagrange  point  L3  is  determined  as 
y  =  .9999956277060894  in  nondimensional  units.  Without  solar 

3 

radiation  pressure,  the  location  of  L3  is  calculated  as 
?3  =  .9999982264196291.  Hence,  the  location  of  L3  in  the  CR3BP  moves 
approximately  388.762012  kilometers  closer  to  the  Sun  when  solar 
radiation  pressure  forces  are  included  in  the  model. 

The  addition  of  the  solar  radiation  pressure  force  then  causes  the 
location  of  the  Lagrange  point  to  move  closer  to  the  Sun  for  all  three 
collinear  libration  points.  It  would  also  be  logical  to  question  how 
solar  radiation  pressure  forces  affect  the  locations  of  the  remaining 
two  Lagrange  points.  The  existence  and  locations  of  the  triangular 
Lagrange  points  when  solar  radiation  pressure  forces  are  included  may 
be  valuable,  but  such  an  investigation  is  not  part  of  this  effort.  The 
focus  here  now  turns  to  the  locations  of  the  Lagrange  points  in  the 
ER3BP,  which  is  the  dynamic  model  used  for  this  work. 
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2.  The  ER3BP 


Five  libration  points  also  exist  in  the  ER3BR,  but  they  arc  not 
stationary  relative  to  the  rotating  frame;  rather,  the  col  linear  points 
pulsate  along  the  xr  axis,  and  the  triangular  points  pulsate  relative 
to  both  the  x  and  the  y  axes  as  the  distance  between  the  primaries 

R  R 

varies  with  time.  The  equilibrium  solutions  can  be  located  by  using 
equations  (2-26)  through  (2-31)  to  find  ratios  of  certain  distances 
that  are,  in  fact,  constant  in  the  problem.  The  collinear  libration 
points  in  the  ER3BP  can  be  found  by  assuming  x  *  0,  x  *  0,  and  y  =  y  = 
z  =  z  =  y  =  z  =  0. 

With  these  restrictions  on  the  coordinates  of  the  collinear 
points,  equation  (2-27)  reduces  to 


20x  +  0x  =  O. 


(2-48) 


Also,  equations  (2-26)  and  (2-29)  combine  to  produce 


(l-/i)(x+Rp)  p(x-R+pR)  ks(x+Rp) 


+  02x, 


(2-49) 


where 


d  =  I  x  +  Rp  I , 


(2-50) 


r  =  |  x  —  R  +  pR  | , 


(2-51 


R  =  1  -  e  cos(E) . 


(2-52) 


Equations  (2-48)  and  (2-49)  can  be  combined  to  provide  one 
equation  written  in  terms  of  the  coordinate,  x,  of  a  particular 
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collinear  Lagrange  point.  Solving  for  x  in  equation  (2-48)  results  in 


(2-53) 


Expressions  for  0  and  0  must  be  substituted  into  equation  (2-53). 
Using  equations  (2-32),  (2-33),  and  (2-52),  equation  (2-53)  simplifies 
to 


e  sin(E)  x  _  e  sin(E)  x 
(1-e  cos (E) )2  R2 


(2-54) 


The  derivative  of  equation  (2-53)  with  respect  to  nondimensional  time 
produces  the  equation 


x 


e  cos(E)  -  e2  x 


(2-55) 


for  which  the  scaled  version  of  equation  (2-i8),  that  is 


E  =  — ,  (2-56) 

R 

is  required. 

Using  equations  (2-32)  and  (2-55)  in  equation  (2-49),  the  term 

..  . 

(x  -  0'x)  reduces  to 
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(X  -  e2x)  = 


(2-57) 


x 


R 


Finally,  using  equation  (2-57)  in  equation  (2-49),  an  expression  for 
the  x  coordinate  associated  with  any  of  the  three  coll  inear  Lagrange 
points  in  the  ER3BP  can  be  derived: 


— -  d3r3+  (p-l)(-^-  +  p)r3-  p(-2 — (2-p))d3+  ks(-j^-  +  lt)r3=  0-  (2-58) 


Expressions  for  x  as  a  function  of  y  (corresponding  to  each  collinear 
libration  point)  can  then  be  substituted  into  (2-58)  to  find  numerical 
values  for  y  ,  y  ,  and  y^.  For  the  ER3BP,  these  expressions  are 


x  =  Li  =  R(l-p-2ri), 
x  =  Lz  =  R  (  1  -  p  +  7T2) , 
x  =  L3  =  R  (  -  n  -  r3  )• 


(2-59) 

(2-60) 

(2-61) 


The  relative  locations  of  the  Lagrange  points  in  the  ER3BP  are  depicted 
in  Figure  2-3. 

When  equation  (2-59),  (2-60),  or  (2-61)  is  substituted  into 
equation  (2-58),  equations  (2-43),  (2-44),  or  (2-45),  respectively, 
result.  Therefore,  the  relative  positions  of  the  collinear  libration 
points  can  be  calculated  from  the  same  three  quintics  that  were  solved 
in  finding  the  locations  of  the  collinear  points  in  the  CR3BP.  The 
locations  of  the  three  points  can  then  be  calculated  at  any  time  by 
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substituting  the  appropriate  numerical  values  of  y  ,  z2>  and  y3  into 
equations  (2-59),  (2-60),  and  (2-61),  respectively.  Although  ^  are 
constants,  the  locations  for  the  collinear  points  on  the  rotating  x 
axis  vary  with  time  because  R(t)  is  time-varying,  ihe  value  of  e  used 
in  these  calculations  is  .01671173051. 

For  this  work,  equation  (2-59)  is  used  to  determine  the  position 
of  relative  to  the  barycenter,  B;  it  is  also  used  to  transform  the 
numerically  integrated  data  from  rotating  coordinates  with  origin  at 
the  libration  point  to  rotating  coordinates  centered  at  the  barycenter. 


E.  State  Transition  Matrix 

For  determination  of  an  acceptable  nominal  trajectory,  as  well  as 
orbit  determination  error  analysis  investigations  and  station-keeping 
studies,  it  is  required  that  the  state  transition  matrix  be  available 
at  predetermined  but  varying  time  intervals  along  the  nominal  path. 
The  transition  matrix  is  derived  in  connection  with  a  linearizing 
analysis. 

The  equations  of  motion  for  the  infinitesimal  mass  in  the  ER3BP 
can  be  linearized  about  a  reference  trajectory  (nominal  path)  that  is  a 
solution  of  the  differential  equations.  The  states,  three  position  and 
three  velocity,  and  the  state  vector  x  are  defined  as 


x  =  x,  x  =  y ,  x  =  z 
i  2  3  3 


x 

4 


(2-62) 


and 


x  =  (x 


(2-63) 


The  reference  trajectory  is  defined  as  x  .  Therefore,  using  a 

REF 

Taylor’s  series  approach,  the  expansion  about  the  reference  path  is 
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written  in  the  form  of  the  fi~st-order  variational  equation 


(  x)  =  x  =  A ( t )  x 

dt 


(2-64) 


where  x  =  x  -  x  is  understood  to  be  the  vector  of  residuals  relative 

REF 

to  the  nominal  solution,  and  the  matrix  A(t)  contains  the  first-order 
terms  in  the  Taylor’s  series  expansion  of  the  equations  of  motion  about 
the  nominal  or  reference  solution  of  interest. 

Using  equations  (2-26)  through  (2-31),  A ( t )  can  be  expressed  as 


A  ( t )  = 


Urr+Gft  20ft 


where  all  four  submatrices  are  dimension  3x3  and 


Uxx  Uxy  Uxz 
Urr  =  Uyx  Uyy  Uyz 


I  Uzx  Uzy  Uzz  I 


(2-65) 


(2-66) 


0  1  0 

n  =  -i  o  o  . 

0  0  0 


The  entries  of  the  symmetric  matrix  Urr  are  given  by 


ks+ji-1  n  3(  1-p-ks)  (x+Rp)2  3p(x-R(  1  -p) ) 

Uxx= - +  62  +  -  +  - ,  (-'  <■’  ■’ 

-.3  3  ,s  5 

dr  d  r 


ks+p-1  n  3(l-p-ks)y2  3py2 

(Jyy= - -  +  02  +  -  +  - , 

3  3  S  S 

dr  dr 


(2-66 } 


?  2 

ks+p-1  p  3(l-p-ks)z'  3pz 


Uzz— 


,3  3 

d  r 


s  ’ 


(2-69) 


3(l-p-ks) (x+Rp)y  3p(x-R( 1-p) )y 


Uxy= 


(2-70 • 


3(l-p-ks) (x+Rp)z  3p(x-R(l-p) )z 


Uxz= 


(2-71 


3(l-p-ks)yz  3pyz 


Uvz~ 


(2-72) 


The  matrix  A ( t )  can  then  be  evaluated  along  the  reference  trajectory 

The  vector  differential  equation  (2-64)  governing  t  ho  n!-i‘e 
variations  from  the  nominal  path  has  a  solution  of  the  form 
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(2-73) 


x(t)  =  <Ht,  tQ)  x(tQ) 


where  4> ( t ,  tQ )  is  the  state  transition  matrix  at  Lime  "t"  relative  to 

time  "t  The  matrix  4>,  then,  represents  the  sensitivities  of  the 
o 

states  at  time  "t"  to  small  changes  in  the  initial  conditions.  It  is 
determined  by  numerically  integrating  the  matrix  differential  equation 

~  4>( t,  t  )  =  <J>{  t ,  t  )  =  A  ( t )  4>(  t ,  t  ),  (2-74) 

dt  o  o  0 


with  initial  conditions  <J>(t  ,t  )  =  I,  the  6x6  identity  matrix.  Thus, 

o  o 

the  nonlinear  equations  of  motion  in  (2-26)  through  (2-28)  and  the 
matrix  equation  (2-74)  combine  to  result  in  42  first-order  differential 
equations  that  can  be  simultaneously  integrated  numerically  to 
determine  the  state  vector  and  its  associated  transition  matrix  at  any 
instant  of  time.  The  reference  trajectories  that  are  of  interest  in 
this  research  are  generated  by  a  numerical  integration  method  that  uses 
a  differential  corrections  process  developed  by  Howell  and 
Pernicka. 112,26-301  The  orbits  include  solar  radiation  pressure  forces 

(31 J 

as  formulated  by  Bell  specifically  for  an  orbit  associated  with  the 

interior  Lagrange  point  in  the  Sun-Earth  system.  The  numerical 

integration  routines  used  in  this  work  are  fourth-  and  fifth-order 

(32J 

Runge-Kutta  formulas  available  in  the  386-Matlab  software  package. 


F.  Bounded  Orbits  Near  Libration  Points 

The  computation  of  bounded  periodic  and  quasi -periodic  orbits  in 
the  vicinity  of  libration  points  has  been  of  increasing  interest  during 
the  past  100  years.  The  accomplishments  of  Poincare  (1854-1912) 
include  a  significant  body  of  work  on  fami I i'-s  of  bounded  periodic 
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orbits  in  the  planar  CR3BP,  the  planar  ER3BP,  and  the  three-dimensiono 
ER3BP.  During  the  time  per  od  1913  to  1939,  Professor  Elis  Stromgre-. 
and  others  at  the  Copenhagen  Observatory  computed  a  family  of  bounde  ! 
periodic  orbits  for  p=.50  near  a  col  linear  libration  point  in  Lh  - 
planar  CR3BP.  141  Other  researchers  have  constructed  two-dimensional 
and  three-dimensional  bounded  periodic  and  quasi-pei iodic  orbits  near 
Lagrange  points;  some  have  investigated  stable  periodic  orbits  in  '  1>. 

-  ...  .  ..  ..  .  .  11,2,5-8,12-22,25-31) 

vicinity  of  the  coll  inear  libration  points. 

This  section  first  discusses  the  stability  of  the  libration  points 
in  the  CR3BP  and  then  the  ER3BP.  The  construction  of  bounded  orbits 
near  the  collinear  Lagrange  points  is  then  summarized.  Finally,  the 
specific  reference  trajectories  used  in  the  orbit  determination  error 
analysis  and  station-keeping  studies  in  this  work  are  introduced. 


1.  Stability  of  the  Libration  Points  in  the  CR3BP 

The  accomplishments  of  those  researchers  who  have  constructed 
bounded  orbits  near  collinear  libration  points  are  particularly 
significant  because  the  collinear  points  are  considered  "unstable" 
points  of  equilibrium  but  with  (only)  one  mode  producing  positive 
exponential  growth.  Bounded  motion  in  their  vicinity,  therefore,  is 
determined  by  deliberately  not  exciting  the  unstable  mode.  A  second 
mode  produces  negative  exponential  orbital  decay  and  is  also 
deliberately  not  excited.  In  the  CR3BP,  the  remaining  four  eigenvalues 
are  purely  imaginary.  The  existence  of  initial  conditions  that  result 
in  only  trigonometric  (sinusoidal)  functions  as  solutions  means  that 
the  collinear  libration  points,  while  unstable,  possess  conditional 
stability  (with  proper  choice  of  initial  conditions)  in  the  linear 

(9) 

sense. 
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The  instability  of  the  collinear  libration  points  can  be 
demonstrated  by  looking  at  a  specific  example  formulated  in  the  CR3BP. 
Earlier  in  this  chapter,  the  equations  of  motion  in  the  ER3BP  were 
linearized  about  a  nominal  path  to  give  the  first-order  variational 
equation  (2-64)  that  includes  the  Jacobian  matrix  A(t).  The  entries  of 
the  matrix  A { t )  are  then  defined  in  equations  (2-65)  through  (2-72). 
These  equations  can  be  simplified  by  assuming  the  CR3BP  is  valid; 
hence,  0=0  and  R  =  8  =  1.  Also,  if  the  stability  of  the  collinear 
points  is  to  be  determined,  the  reference  solution  can  now  be  the 
(stationary)  libration  point  of  interest  in  the  CR3BP.  Under  the 
assumptions  of  the  CR3BP  and  the  libration  point  reference  solution, 
the  Jacobian  matrix  Aft)  is  constant  and  can  be  denoted  as  the  constant 
matrix  A;  stability  measures  such  as  eigenvalue  analysis  can  then  be 
accompl ished. 

For  this  analysis,  the  vector  x  represents  relatively  small 

variations  from  a  particular  libration  point  of  interest  (L  ).  The 

i 

linear  variational  equation  consequently  may  be  written  as 


x  =  A  x 


(2-74) 


where  x  =  (a  (3  y  a  (Z  ?]  ,  and  A  is  the  constant  Jacobian  matrix.  In 
order  to  compute  the  constant  matrix  A,  nondimensional  collinear 
libration  point  coordinates  (x  ,0,0)  for  the  CR3BP  are  substituted 
into  the  matrix  A.  (The  matrix  Aft)  is  given  by  equations  (2-65) 
through  (2-72),  and,  under  the  assumptions  that  the  CR3BP  is  valid  and 
the  given  libration  point  is  the  nominal  solution,  it  reduces  to  the 
constant  matrix  A.)  The  matrix  A  will  obviously  be  numerically 
different  for  each  of  the  five  libration  points  of  interest. 


The  following  .  \ree  linear  equations  result: 


a  =  2  0  +  (2B+1 )  a 


(2-75) 


0  =  (2-B)0  -2a 


(2-76) 


i  =  (-B)r 


where 


1-ji+ks  p 

[XLl+,i]3  [XLi+^"1]3 


(2-77) 


(2-78) 


and  x  represents  the  nondimens ional  coordinate  of  the  coll inear 
Li 

libration  point  of  interest.  The  variational  equation  for  the 
out-of-plane  component  y  decouples  from  the  in-plane  equations  and 
represents  pure  oscillatory  motion.  The  equations  in  terms  of  a  and  (3 
produce  the  characteristic  equation 

A4  +  (2-B)A2  +  (1+B-2B2)  =  0  (2-79) 

2 

with  roots  that  depend  on  the  sign  of  (2-B+2B  ).  For  the  three 
collinear  Lagrange  points,  the  sign  of  this  term  is  negative  for  all 
values  of  the  mass  ratio  such  that  Osps.5.  Consequently,  the 
characteristic  equation  (2-79)  has  roots  that  include  a  purely 
imaginary  pair  and  a  real  pair,  opposite  in  sign.  The  solution  of  the 
linear  variational  equations  (2-75),  (2-76),  and  (2-77)  for  any  of  the 
collinear  points  then  shows  unbounded  motion  (exponential  growth)  in 
response  to  arbitrary  initial  conditions.  The  benefits  of  selecting 
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specific  initial  conditions  to  avoid  exciting  the  unbounded  exponential 
modes  will  be  discussed  later  in  this  section.  But  first,  a  short 
paragraph  discussing  stability  of  the  triangular  Lagrange  points  in  the 
CR3BP  is  included  because  the  discoveries  of  natural  satellites 
orbiting  near  triangular  libration  points  has  helped  spur  the  current 
substantial  interest  in  libration  point  spacecraft  orbits. 

The  triangular  libration  points  are  marginally  stable  in  the 
linear  sense  for  a  specific  range  of  primary  mass  ratio  in  the  CR3BP. 
Purely  imaginary  roots  in  two  conjugate  pairs  are  obtained  for  ps.0385, 
which  is  given  here  to  four  decimal  places  and  is  sometimes  referred  to 
as  Routh’ s  value.141  The  mass  ratios  (listed  here  to  three  decimal 
places),  for  example,  in  the  three-body  systems  of  the  Earth-Moon 
(p  =  1.216  x  10  6),  Sun-Ear th+Moon  (p  =  3.022  x  10  6)  and  Sun-Jupiter 
(p  =  9.485  x  10~4)  all  satisfy  the  mass  ratio  requirement  for  marginal 
stability  of  the  triangular  points  in  the  linearized  model.  Natural 
satellites,  such  as  the  Trojan  asteroids  or  a  moon  of  Saturn,  occupy 
linearly  stable  orbits  near  triangular  libration  points  in  their 
respective  systems. 


2.  Stability  of  the  Libration  Points  in  the  ER3BP 

Several  researchers  have  analyzed  the  stability  of  the  libration 
points  in  the  elliptic  problem,  where  both  the  mass  ratio,  p,  and  the 
primary  orbit  eccentricity,  e,  influence  stability. t4'9’33-351  The 
instability  of  the  collinear  libration  points  as  determined  in  the 
circular  problem  for  all  the  values  of  mass  parameter  persists  for  the 
elliptic  problem;  an  analysis  of  the  collinear  points  shows  instability 
for  any  combination  of  the  values  of  both  p  and  e. 

The  results  of  a  linearized  stability  analysis  regarding  the 
effects  of  eccentricity  and  mass  ratio  on  the  linear  stability  of  the 

[34  ] 

triangular  points  have  been  published  by  Danby  and  then  later  by 
[35] 

Bennett  .  Both  Danby  and  Bennett  have  numerically  generated  graphic 
depictions  of  the  linear  stability  region  in  the  p-e  plane.  For  the 
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eccentricity  in  the  Sun-Ear th+Moon  ER3BP,  the  value  of  ; i  which  ensui <  , 
linear  stability  is  only  slightly  less  than  Routh’s  value  (decreased  i 
approximately  one  percent).  An  interesting  aspect  of  the  p-e  stabilit  j 
region  is  that  a  range  of  values  of  p  greater  than  Routh’s  value  al.o 
defines  a  region  of  linear  stability  for  a  specific  range  of  values  cl 
e  less  than  .3143. 


3.  Construction  of  Bounded  Collinear  Libration  Point  Orbits 

The  initial  goal  in  the  process  of  generating  bounded  orbits  near 
a  collinear  (unstable)  libration  point  is  to  avoid  exciting  the 
unstable  mode  associated  with  the  linearized  motion.  The  meteoric  dust 
particles  that  may  be  orbiting  near  Lagrange  point  in  the  Sun-Earth 
system  could  only  linger  near  that  point  if  they  arrive  with  the 
"correct"  initial  position  and  velocity  states  relative  to  L^.  The 
"correct"  initial  conditions  will  only  (primarily)  excite  the  stable 
modes  associated  with  the  linearized  motion  and  not  (or  minimally) 
excite  the  unstable  mode.  The  degree  to  which  the  unstable  mode  is 
excited  will  determine  the  length  of  time  that  the  dust  particles 
linger  near  L ,  .  The  general  solution  to  the  harmonic  equation  (2-77) 
is  written 

T'(t)  =  B  cos(ut)  +  B  sinCi/t)  (2-81) 

1  2 

where  Bt  and  are  constants  of  integration,  with  values  determined  b/ 

initial  conditions,  and  v  =  (B)  .  The  coupled  in-plane  equations 

have  eigenvalues  of  the  form  iA,  -iA,  77,  -tj  for  A>0  and  jj>0  where 
1  /2 

i  =  (-1)  .  They,  therefore,  yield  solutions 


a(t)  =  A^e 


4t 


A  e 
2 


-7)t 


A  cos (At) 
3 


A  sin(Atl 

4 


(2-82) 
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(2-83) 


(3(  t)  =  -C  A  e 


7}t 


C  A  e 
1  2 


“T)t 


C  A  cos(At) 

2  4 


-  C  A  sin(At) 
2  3 


where 


C 

l 


—  (1  +  2B  -  u2),  C  =  _L  (1  +  2B  +  A2) 
2r)  2  2A 


(2-84) 


and  A  ,  A  ,  A  ,  A  are  constant  coefficients,  dependent  on  the  given 
initial  conditions.  For  bounded  motion,  the  coefficient  A^  must  be 
zero.  For  both  periodic  and  quasi-periodic  motion,  A  must  also  be 
zero.  These  two  restrictions  eventually  lead  to  two  equations  that  the 
initial  conditions  must  satisfy  in  order  to  produce  a  linear  solution 
that  is  bounded.  The  initial  conditions  are  denoted  aQ,  /3q,  ao>  /3q. 

The  constant  coefficients  (A  ,A  ,A  ,A  )  can  then  be  treated  as 

1  2  3  4 

unknowns  and  expressed  in  terms  of  the  initial  conditions.  Both  A  and 

a  -11  u  .  r  [4 ) 

A  will  be  zero  if 
2 

A  =  F  a  +  F  P  +  F  a  -  .5/3  =  0  (2-85) 

1  1  0  2  0  3  0  o 


A  =  F  a 
2  l  o 


F  13 
2  0 


.5  (i  =  0 
o 


(2-86) 


where 


-A 

F  =  - - -  ;  F  =  - - -  ;  and  F  =  .5  C  .  (2-87) 

1  2 (AC  — t}C  )  2  2(tjC  +AC  )  3  2 

2  1  2  1 


When  equations  (2-85)  and  (2-86)  are  satisfied,  the  sinusoidal 
solutions  to  the  linear  equations  (2-75),  (2-76),  and  (2-77)  can  be 
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written  in  more  compact  form  (using  trigonometric  identities)  as 


a(t)  =  -A  cos(At  +  <p) ,  (2-89) 

X 

(3(t)  =  C  A  sin(Xt  +  <(>),  (2-90) 

2  x 

y(t)  =  A  sin(i/t  +  4>) ,  (2-91) 

2 

A  =  (A2  +  A 2  )1/2,  (2-92) 

x  3  4 

A  =  (B2  +  B2  )1/2,  (2-93) 

2  1  2 

tan($)  =  -  (A  /  A  ),  (2-94) 

4  3 

tan(^r)  =  (B  /  B  ) .  (2-95) 

12 


The  linear  solution  can  be  very  useful  in  developing  a  higher- 

order  approximation  to  the  actual  nonlinear  trajectory.  The  motion 

represented  by  the  linearized  solution  as  it  appears  above  will  be 

quasi-periodic  since  A  and  v  are,  in  general,  not  equal.  They  are, 

however,  of  the  same  order  of  magnitude  for  the  problem  of  interest. 

Periodic  motion  can  be  constructed  by  modifying  equation  (2-77)  to 

(17) 

force  a  solution  with  equal  in-plane  and  out-of-plane  frequencies: 
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a  =  2  0  +  (2B+1 )  a 


(2-96) 


0  =  (2-B)0  -2a  (2-97) 

i  +  \\  =  0.  (2-98) 


Equation  (2-98)  is  derived  by  adding  (A2  -  B)y  to  the  left-hand 

side  of  equation  (2-77);  the  addition  of  this  term  is  compensated  for 

by  using  the  higher-order  terms  in  the  analytic  approximation. 

Specifically,  periodic  linear  motion  can  be  obtained  if  the  in-plane 

(A^)  and  out-of-plane  (A^)  amplitudes  of  the  motion  are  related  in  such 

a  way  that  allows  nonlinear  contributions  (higher-order  terms  in  the 

analytic  approximation)  to  the  system  to  produce  eigenfrequencies  that 

[17] 

are,  in  fact,  equal.  The  resulting  linearized  periodic  equations 

would  be  of  the  form 


a(t)  =  -A  cos(At  +  0) 

X 

(2-99) 

0(t)  =  C  A  sin(At  +  <p) 

2  x 

(2-100) 

y(t)  =  A  sin(At  +  ip) . 

(2-101) 

The  third-order  analytic  representation  is  used  in  this  work  to 
provide  the  initial  model  for  the  trajectories.  The  method  of 
successive  approximations,  using  the  linearized  solution  as  the  first 
approximation  to  the  nonlinear  orbital  path,  and  the  method  of  dual 

f IS  |s  25 j 

time  scales  are  used  to  derive  the  third-order  result.  ’  ’  The 

method  of  successive  approximations  is  used  to  generate  an  asymptotic 
series  in  an  appropriately  small  parameter.  (The  square  root  of  the 
eccentricity  of  the  primary  orbit,  that  is  the  orbit  of  the  Earth-Moon 


52 


barycenter  about  the  Sun,  is  the  small  parameter  used  here. )  Th 
method  of  dual  time  scales  is  used  to  convert  the  system  of  ordinal  . 
differential  equations  into  a  system  of  partial  differential  equations 
In  general,  the  method  of  multiple  scales  allows  the  various  nonlinc.f 
resonance  phenomena  to  be  included  in  the  approximate  analytic  solulk n 
and  provides  a  method  to  remove  secular  terms.  (Here,  "secular"  refers 
to  terms  that  include  the  time  variable  and  is  derived  from  the  French 
"siecle"  meaning  century. ) 

(25) 

The  analytic  solution  of  Richardson  and  Cary  for  the 

Sun-Earth+Moon  ER3BP  has  been  derived  to  fourth  order,  but  the  third- 

order  approximation  is  found  to  be  sufficient  for  this 

research.  ’  A  numerical  integration  algorithm,  using  a 

differential  corrections  procedure  that  is  designed  to  adjust  the  first 

guess  as  obtained  from  the  analytic  approximation,  can  then  be  used  to 

numerically  generate  the  orbit  of  interest.  A  method  developed  by 

Howell  and  Pernicka  ’  is  used  here  to  generate  the  orbital 

paths.  Their  method  initially  employs  the  approximate  analytic 

solution  to  compute  target  points.  A  two-level  (position  matching  then 

velocity  matching),  multi-step  differential  corrections  algorithm  is 

used  to  construct  a  numerically  integrated,  bounded  trajectory  that  is 

continuous  in  position  and  velocity.  A  solar  radiation  pressure  model 
[31 J 

developed  by  Bell  is  also  incorporated  in  the  numerical  integration 
procedure. 

The  method  of  Howell  and  Pernicka,  including  solar  radiation 
pressure,  uses  an  initial  analytic  guess  that  represents  a  halo  orbit 
or,  alternatively,  a  considerably  smaller  Lissajous  path.  The 
higher-order  terms  tend  to  slightly  alter  the  first-order  periodic  or 
quasi-periodic  path.  Consequently,  the  initial  target  path  for  a  halo 
orbit  will  generally  not  be  precisely  periodic.  The  two-level, 
multi-step  differential  corrections  procedure  then  adjusts  the  initial 
analytic  target  orbit  and,  therefore,  will  compute  a  halo-type  orbit 
that  is  nearly  (but  not  exactly)  periodic. 
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4.  The  Reference  Paths  Generated  for  This  Work 


An  original  goal  of  this  research  was  to  compare  orbit 

determination  error  analysis  results  and  station-keeping  costs  for 

Lissajous  and  halo  orbits.  Precisely  periodic  halo  orbits  exist  in  the 

CR3BP .  They  also  exist  in  the  ER3BP,  but,  in  the  ER3BP ,  they  are 

multiple  revolution  trajectories  with  periods  much  longer  than  those  of 

interest  here.  Nearly  periodic  orbits  are  more  practical  in  the  ER3BP 

and  are  much  more  likely  to  be  used  in  mission  planning;  therefore,  the 

goal  here  should  be  slightly  modified  to  be  the  comparison  of  Lissajous 

and  "halo-type"  orbits.  The  general  shapes  of  the  three-dimensional 

halo-type  and  Lissajous  orbits  can  be  seen  by  plotting  three 

orthographic  views  of  each  orbit,  using  the  tabular  data  from  the 

numerical  integration  routine.  Figure  2-4  depicts  three  orthographic 

views  of  point  plots  for  the  Lissajous  orbit  used  in  this  research. 

Figure  2-5  contains  three  orthographic  views  (on  a  slightly  different 

scale)  of  the  considerably  larger  halo-type  orbit  to  be  used  later  in 

this  work.  (Note  that,  in  general,  the  amplitude  ratio  for  Lissajous 

trajectories  is  arbitrary.  In  halo  orbits,  however,  constraining  the 

amplitude  ratio  results  in  equalized  frequencies  for  in-plane  and 

out-of-plane  motion. )  The  orbits  are  depicted  in  the  rotating 

reference  frame  centered  at  L  . 

l 

Both  orbits  are  clearly  not  periodic;  a  Lissajous  orbit  is  often 
called  a  quasi-periodic  path,  and  these  two  orbits  could  clearly  be 
termed  quasi-periodic  or  Lissajous  paths.  The  major  difference  between 
the  orbits  is  the  larger  size  of  the  halo-type  orbit;  however,  other 
differences  are  also  present.  The  maximum  x  and  y  excursions  of  the 
halo-type  orbit  are  approximately  four  times  as  large  as  those  of  the 
Lissajous  path.  Furthermore,  the  direction  of  motion  (clockwise  versus 
counterclockwise),  as  viewed  in  the  y-z  orthographic  depiction,  is 
different  for  the  two  orbits  used  here.  The  direction  of  motion  on  the 
halo-type  orbit  is  counterclockwise  in  the  y-z  depiction;  the  direction 
of  motion  is  clockwise  in  the  y-z  depiction  for  the  Lissajous  path. 
(Both  orbits  include  clockwise  motion  in  the  x-y  depiction.) 
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xios  „  Xl05 


x  in  kilometers  xl04  x  in  kilometers  xlO4 


Indicates  the  direction  of  motion. 


Indicates  the  location  of  the  libration  point. 


Figure  2-4.  Three  Orthographic  Views  of  a  Lissajous  Orbit. 
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y  in  kilometers  xl()5 


indicates  direction  of  motion. 

.  indicates  the  libration  point  location. 


Figure  2-5.  Three  Orthographic  Views  of  a  Halo-Type  Orbit. 
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The  two  orbits  can  also  be  differentiated  in  terms  of  In¬ 

direction  of  the  maximum  z  excursion  in  the  x-z  depiction.  If  the 
maximum  z  excursion  is  in  the  positive  z  direction,  the  orbit  can  be 
termed  a  member  of  a  "northern  family"  of  orbits.  When  the  maximum 
excursion  of  the  orbit  is  in  the  negative  z  direction,  the  orbit  is 

termed  a  member  of  a  "southern  family”  of  orbits.  In  the  x-z 

orthographic  depiction,  the  smaller  (Lissajous)  path  can  be  seen  to  be 
a  member  of  a  northern  family  of  orbits,  while  the  halo-type  orbit  is  a 
member  of  a  southern  family  of  orbits. 

The  differences  in  the  paths,  outlined  in  the  preceeding 

paragraphs,  are  the  consequences  of  the  initial  conditions  selected  for 
the  numerical  integration  of  the  trajectory.  The  initial  conditions 
used  for  computation  of  the  orbit  near  libration  point  L  determine  the 
related  characteristics  of  direction  of  motion  in  the  y-z  plane  and  the 
family  of  orbits  (northern  or  southern)  to  which  it  belongs. 


CONCLUSION 


The  restricted  three-body  problem  is  an  important  and  interesting 
area  of  research.  The  Sun-Earth+Moon  ER3BP,  in  particular  is  currently 
an  area  of  vital  research  attention.  The  addition  of  solar  radiation 
pressure  in  the  Sun-Earth+Moon  ER3BP  alters  the  locations  and  the 
conditions  for  linear  stability  of  the  Lagrange  points.  Including 
these  solar  effects  improves  the  accuracy  of  the  dynamic  model,  but 
there  are  also  additional,  related  benefits.  Simulations  used  to 
calculate  expected  tracking  errors  and  the  cost  of  maintaining  the 
spacecraft  near  the  unstable  orbits  can  now  also  include  uncertainty  in 
this  solar  radiation  pressure  force.  The  addition  of  this  real 
uncertainty  will,  in  turn,  improve  the  accuracy  of  further  studies 
conducted  in  relation  to  future  spacecraft  missions. 
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